﻿using System;
using System.Collections.Generic;
using Allegro.Integrator;


namespace Cosmologica.Integrator.RungeKutta
{
    /*
    This is the Butcher tableau
    Define f = 1398169080000
    db = b-bh
    c		a*f
    ===========
    0
    1/6		233028180000
    4/15	74569017600	    298276070400
    2/3		1165140900000	-3728450880000	3495422700000
    5/6		-3604654659375	12816549900000	-9284716546875	1237962206250
    1		3355605792000	-11185352640000	9172628850000	-427218330000	482505408000
    1/15	-770204740536	2311639545600	-1322092233000	-453006781920	326875481856	0
    1		2845924389000	-9754668000000	7897110375000	-192082660000	400298976000	0		    201586000000

    b*f		104862681000	0		        545186250000	446637345000	188806464000	0		    15076875000	    97599465000

    db*f	8738556750	    0		        9735468750	    -9709507500	    8582112000	    95329710000	-15076875000	-97599465000
     */

    /// <summary>
    /// This is a stage 8, order (5,6) Runge-Kutta pair published by Jim Verner.
    /// This is the RK pair implemented by the dverk subroutine of CMBFAST.
    /// 
    /// </summary>
    public class RKTableau_56_CMBFAST : ButcherTableau
    {
        public RKTableau_56_CMBFAST()
            : base(8)
        {
            double f = 1398169080000.0;

            // Initialize the vectors
            c[0] = 0.0;
            c[1] = 1.0 / 6.0;
            c[2] = 4.0 / 15.0;
            c[3] = 2.0 / 3.0;
            c[4] = 5.0 / 6.0;
            c[5] = 1.0;
            c[6] = 1.0 / 15.0;
            c[7] = 1.0;

            a[1, 0] = 233028180000.0 / f;

            a[2, 0] = 74569017600.0 / f;
            a[2, 1] = 298276070400.0 / f;

            a[3, 0] = 1165140900000.0 / f;
            a[3, 1] = -3728450880000.0 / f;
            a[3, 2] = 3495422700000.0 / f;

            a[4, 0] = -3604654659375.0 / f;
            a[4, 1] = 12816549900000.0 / f;
            a[4, 2] = -9284716546875.0 / f;
            a[4, 3] = 1237962206250.0 / f;

            a[5, 0] = 3355605792000.0 / f;
            a[5, 1] = -11185352640000.0 / f;
            a[5, 2] = 9172628850000.0 / f;
            a[5, 3] = -427218330000.0 / f;
            a[5, 4] = 482505408000.0 / f;

            a[6, 0] = -770204740536.0 / f;
            a[6, 1] = 2311639545600.0 / f;
            a[6, 2] = -1322092233000.0 / f;
            a[6, 3] = -453006781920.0 / f;
            a[6, 4] = 326875481856.0 / f;
            a[6, 5] = 0.0;

            a[7, 0] = 2845924389000.0 / f;
            a[7, 1] = -9754668000000.0 / f;
            a[7, 2] = 7897110375000.0 / f;
            a[7, 3] = -192082660000.0 / f;
            a[7, 4] = 400298976000.0 / f;
            a[7, 5] = 0.0;
            a[7, 6] = 201586000000.0 / f;

            // High order weights
            b[0] = 104862681000.0 / f;
            b[1] = 0.0;
            b[2] = 545186250000.0 / f;
            b[3] = 446637345000.0 / f;
            b[4] = 188806464000.0 / f;
            b[5] = 0.0;
            b[6] = 15076875000.0 / f;
            b[7] = 97599465000.0 / f;

            // Low order weights
            bh[0] = 96124124250.0 / f;
            bh[1] = 0.0;
            bh[2] = 535450781250.0 / f;
            bh[3] = 456346852500.0 / f;
            bh[4] = 180224352000.0 / f;
            bh[5] = -95329710000.0 / f;
            bh[6] = 30153750000.0 / f;
            bh[7] = 195198930000.0 / f;
        }
    }
}
